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Abstract 

We develop an agent-based model of the motion and pattern formation of vesicles. These 
intracellular particles can be found in four different modes of (undirected and directed) 
motion and can fuse with other vesicles. While the size of vesicles follows a log-normal 
distribution that changes over time due to fusion processes, their spatial distribution gives 
rise to distinct patterns. Their occurrence depends on the concentration of proteins which are 
synthesized based on the transcriptional activities of some genes. Hence, differences in these 
spatio-temporal vesicle patterns allow indirect conclusions about the (unknown) impact of 
these genes. 

By means of agent-based computer simulations we are able to reproduce such patterns on 
real temporal and spatial scales. Our modeling approach is based on Brownian agents with an 
internal degree of freedom, 6, that represents the different modes of motion. Conditions inside 
the cell are modeled by an effective potential that differs for agents dependent on their value 
9. Agent's motion in this effective potential is modeled by an overdampted Langevin equation, 
changes of 9 are modeled as stochastic transitions with values obtained from experiments, 
and fusion events are modeled as space-dependent stochastic transitions. Our results for the 
spatio-temporal vesicle patterns can be used for a statistical comparison with experiments. 
We also derive hypotheses of how the silencing of some genes may affect the intracellular 
transport, and point to generalizations of the model. 

1 Introduction 

Agent-based modeling has proven to be a versatile tool to simulate processes of structure for- 
mation bottom up. By assuming features and interaction rules of agents on the "microscopic" 
level, one is able to observe the emergent systems properties on the macroscopic level. This is 
of particular importance in those areas where the systems dynamics can hardly be captured top 
down, i.e. in living systems, including biological, social or economic systems. 

But the advantage of agent-based models in freely defining agent properties and interactions soon 
turns out to be a pitfall, because this way arbitrary patterns can be generated and it is difficult 
to choose the right values in a high dimensional parameter space. To minimize these problems, 
there are basically two ways: (i) to closely link the agent's properties to experimentally observed 
data, and (ii) to apply methods that allow to aggregate the agent dynamics, to formally derive 
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the systems dynamics. The latter provides a firm relation between agent's features and systems 
feature's and may reveal also the role of certain (control) parameters. 

The concept of Brownian agents |16) was developed to facilitate the second way. The dynamics of 
agents is described in a stochastic manner, similar to the Langevin approach of Brownian motion. 
This allows to obtain on the macroscopic level a closed-form partial differential equation for the 
density, that for the case of Brownian motion simply describes a diffusion process. The dynamics 
in most real systems, however, is much more complicated. Agents are not simple random walkers, 
they respond to information in their environment, follow chemical gradients, and can at the same 
time also contribute to generating information, chemical gradients etc. Further, agents do not 
behave the same all the time. Instead, they may have different modes of activity each of which 
corresponds to a particular behavior. To cope with these features, Brownian agents are described 
by internal degrees of freedom and their environment is modeled as an adaptive landscape, or 
effective potential, which can be modified by the agents while responding to the information 
provided. Further, transitions between the agent's internal degrees are possible, dependent on 
internal or external conditions. 

On the formal level, the macroscopic dynamics is then no longer described by a diffusion equation, 
but by a quite advanced reaction-diffusion equation with a variable drift term, which is coupled 
to another differential equation describing the dynamics of the adaptive landscape dependent on 
the agent's activity. This allows to tackle the dynamics of systems comprised of many interacting 
agents on two levels: (i) the agent level, where computationally efficient computer simulations can 
be performed, (ii) the system's level, where coupled differential equations may be obtained and 
investigated analytically. The application discussed in this paper is unfortunately complex enough 
to not provide closed form equations for an analytical treatment. Nevertheless, the concept of 
Brownian agents allows us to formally specify the agent dynamics in terms of stochastic equations 
of motion in an adaptive landscape. 

The aim of this paper is twofold: (i) to develop an agent-based model of intracellular transport and 
pattern formation, which is general enough to be applied to various such phenomena involving 
free and directed motion and fusion processes (Section 3), and (ii) to specify this agent-based 
model for the case of vesicle movement and fusion, in close relation to experimental findings 
(Section 4). Importantly, the internal degrees of freedom of agents and transitions between these 
are obtained from experiments. This allows us to observe pattern formation on real time and 
spatial scales (Section 5), the outcome of which can, at least in a statistical manner, be compared 
with real experiments. Hence, applying the concept of Brownian agents to a real problem, i.e. 
the intracellular transport and pattern formation of vesicles, demonstrates both the versatility 
of the concept and its suitability to generate hypotheses about real intracellular processes. 
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2 Vesicle Formation and Vesicle Motion 

In this paper, we are interested in the intracellular transport and pattern formation of vesicles. 
These are quite small intracellular particles (diameter approx. O.l^m) (see |15) . [2] and [6]). 
They are formed at the cell membrane, to contain some extracellular material engulfed by the 
cell membrane. This import of material, called endocytic cargo, may include macromolecules, 
but also viruses, which are all encapsulated in vesicles - a process called endocytosis (see [6], [13j 
and |14)). In this paper, we do not consider endocytosis explicitly, but assume that vesicles are 
formed at the membrane and then released into the interior of the cell at a constant rate (later 
called internalization rate). It is known from experiments that the size distribution of newly 
formed vesicles follows approximately a log-normal distribution (see [T]). 




Figure 1: Two-dimensional representation of a cell with cytoskeletal structure. Adapted from 
[4J. 

Released vesicles can diffuse inside the cell, but they can also be reabsorbed by the membrane 
at a constant rate, later called recycling rate (|3]). Vesicles need to be transported from the cell 
membrane to the endosomal system located in different areas inside the cell - this transport 
process is called membrane trafficking (see j6]). To facilitate this process, in addition to the free 
diffusion^ vesicles can also perform a directed motion along the cytoskeleton, which is an intracel- 
lular structure made of two different kinds of polymer filaments: actin filaments and microtubule 
filaments (see Figure [T]). Whereas actin filaments are randomly distributed, microtubules (MT) 
are directed toward the microtubule organization center (MTOC), which is located close to the 
cell nucleus (see Figure [2]). In order to perform a directed motion along actin or MT filaments, 
vesicles have to bind to motor proteins (kinesines, dyneins and myosins) which basically deter- 
mine the type of motion. The interaction between motor proteins and their related filaments 
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depends on several additional factors, most notably the presence of specific proteins such as Rab 
GTPases, scaffolding proteins and receptor proteins (see [9]). To which extent these proteins are 
present depends on the other hand on the genes of the cell, which have to be transcriptionally 
active in order to synthesize these proteins - a process called gene expression. 

It is exactly this dependency, which motivates our interest in the motion of vesicles. If a given 
gene, for example cdk8^ is transcriptionally active, it allows the synthesis of the protein CDK8, 
which may affect intracellular processes in various, mostly unknown ways. In this paper, we 
precisely ask how such a gene - through the synthesis of the specific protein - affects the motion 
and pattern formation of vesicles. The latter process results from the fact that vesicles can form 
larger ones by fusing with other vesicles. The fusion process relies on energy provided by the cell 
and can only take place on the MT if vesicles are close enough and below a critical size. The 
two simultaneous dynamic processes, namely (free or directed) motion and fusion result in a 
distinct spatio-temporal distribution of vesicles of different sizes - which we want to predict with 
our model. 

The patterns produced by means of our agent-based simulation can then be compared to experi- 
ments which show such vesicle patterns depending on the transcriptional activity of specific genes 
(see [T]). These genes can, for example, be knocked out in RNAi or drug screens, which in turn 
perturb the synthesis of proteins (see [TT]). Of course, such patterns can be only compared in a 
statistical sense, a problem discussed in the Conclusions. However, if we are able to reproduce 
empirical patterns with our model, we argue that the underlying dynamic processes, motion and 
fusion, are covered sufficiently with our modeling assumptions. This does not only hold for the 
assumed interaction between vesicles and MT or other vesicles, it shall also hold for particular 
parameter dependencies, most notably the concentration of specific proteins. Precisely, we want 
to end up with a testable prediction of how these concentrations affect vesicle motion and fusion, 
which shall be confirmed by subsequent experiments (see [1]). Some of the transition rates later 
used in our model specifically depend on the experimental setup, e.g. the endocytic cargo. Here 
we consider the case of Transferrin, an iron-binding protein contained in the vesicles to which flu- 
orescently labelled proteins can be attached, i.e., vesicles can be made visible in the experiment. 
The aforementioned internalization rate and recycling rate are thus taken for Transferrin. 

Eventually, we note that our modeling approach (as every modeling approach) is based on a 
number of simplifications: (i) we neglect the motion of vesicles along actin filaments, because it 
was shown that such processes do not affect the pattern formation (recall that fusion only 
takes place on MT), (ii) we assume that the cytoskeleton is described by the spatial structure 
of MT only (i.e. tha actin filaments are neglected) and that MT are abundant (i.e. there are 
always MT to move on) and do not change in time. This allows to neglect the growth and 
shrinkage of MT when modeling the motion and fusion of vesicles, (iii) We neglect fission, i.e. 
the fragmentation of larger vesicles into vesicles of smaller sizes. 
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3 A Model of Brownian Agents 

Brownian agents Our modeling approach is based on the concept of Brownian agents |16| 
which found many interesting applications in biology [3 |8l [ini [H] but also in modeling social 

systems such as online communities |17) . It allows to formalize the agent dynamics using methods 

(k) 

established in statistical physics. A Brownian agent is described by a set of state variables u\ , 
where the index z = 1, refers to the individual agent while k indicates the different vari- 
ables. These could be either external variables that can be observed from the outside, or internal 
degrees of freedom that can only be indirectly concluded from observable actions. Noteworthy, 
the different (external or internal) state variables can change in the course of time, either due to 
influences of the environment, or due to an internal dynamics. 

In the following, each agent represents an intracellular vesicle which, in accordance with the 
previous description, is able to change its state by spatial mobility, changes of activity, and 
growth processes, as formalized subsequently. 

Spatial mobility For the agent's spatial position rj(t), We assume that changes in the course 
of time can be described by an overdamped Langevin equation of a Brownian particle moving in 
an effective potential |16) : 



The overdamped limit implies that the absolute value of the agent's velocity Vi is approximately 
constant, but the direction may change due to stochastic influences. Further, Vi implicitly depends 
on the agent's mode of activity, Oi{t). 

Eqn. ([T]) assumes that the agent's motion is influenced by two different forces, a deterministic one 
which results from the gradient of the effective potential, and a stochastic one, which is assumed 
to be Gaussian white noise, {ii{t)) = 0, {^i{t)^j{t')) = Sij5{t — t'). The strength of the stochastic 
force D = ksT/'jo determines, in the spatial case, the diffusion coefficient, with 70 being the 
friction coefficient. 

The deterministic part contains two important ingredients: The effective potential h'^(r,t,9) de- 
scribes the conditions inside the cell. The response function a;[^j(t)] depends on the internal state 
of the agent, 9i{t) and determines what component of the effective potential actually influences 
the agent. Both are specified later after we made clear the notion of the internal state 9. 

Changes of activity We assume that the agent's mode of activity 9i{t) can be changed, uj{9'\9) 
being the transition rate from state 9 into any other state 9' . In accordance with the literature 



dt 



Vi[9i{t)] 



amt)] dh%r,t,9) 
70 dr 



+ V2D Ci{t) 
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Free diffusion in the cytosol 

Kinesin-driven transport towards MT plus-ends 

Dynein-driven transport towards MT minus-ends 

Kinesin/Dynein-driven transport on MT with tendency to fusion with other vesicles 
Bound to the cell membrane 



Table 1: Different modes of activities, MT stands for microtubuh 



[1], we distinguish five different modes of activity of a vesicle as shown in Table [T| each of which 
is expressed by a different value of the internal degree of freedom 9. 

While in principle transitions between all states could be assumed, only a subset of them is 
biologically relevant. Table [3] in Section 4.1 will list those together with their respective value, 
i.e. the expected number of transitions per time unit. 



Growth and decay We assume that fusion, i.e. the coalescence of two vesicles with sizes Sj 
and Sj, can be described by a transition rate uj{si-\-j\si, Sj) that depends on the internal states 
6i of the agents, i.e. their ability to fuse, and their effective distance |rj — rj\. As described 
above, fusion is only possible for agents with the internal state 9 = 4 which is expressed by the 
Kronecker delta, 60-^4^. Further, due to the volume exclusion, agents can effectively approach each 
other only up to a distance d (which represents an average spatial extension of vesicles). The 
ability to fuse also depends on the vesicle size because of the energy required for this process. 
Because the available fusion energy is limited, it was observed experimentally that vesicles of 
a size larger than Smax do not fuse. This is considered in the transition rates by an additional 
exponential cut-off term which becomes effectively zero if one of the fusing vesicles reaches the 
maximum size. This leads us to the transition rate for fusion: 

.VI ' J/ d+\ri-rj\ h + e^(^™^-*')] • [1 -h e'^^**'"''^"''^"-'] 

Here uj^ denotes the fusion affinity which depends on the protein concentration, and e = 0.05 is 
chosen as a small number, to increase the cutoff effect. 

Little is known regarding the fission process, i.e. the fragmentation of a vesicle of size si into two 
vesicles of sizes Sj, sj (with si = Si + Sj) . Therefore, we assume that the respective transition rate 
u{si,Sj\si + Sj) is a constant w equal for all possible fragmentation processes, which describes 
spontaneous fragmentation. If w is small compared to other transition rates, fission can be 
neglected in first approximation. 

We note that, because of the fusion process, the total number of vesicles, is no longer constant. 
While a conservation of the total mass M of all vesicles can still be assumed, both the number of 
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vesicles and their size distribution P{Ni, N2, Ns...,t) changes over time. Therefore, we have 

N N 

M = "^3 Ns{t) = const. ; ^ Ns{t) ^ const. (3) 
1=1 1=1 

Effective potential After the above distinction between the different values for the internal 
degree of freedom 6, we can now specify the effective potential which depends on these states. 
h'^{r, t, 6) denotes a scalar potential field that results from the influence of different field compo- 
nents hQ{r,t). Each of these components refers to specific conditions inside the cell. Compared 
to the time scale involved in the motion of the vesicles, some of these conditions can be assumed 
as constant in time, but varying across space. 

With reference to Table [Tj hi{r) describes the influence of the cell membrane on the motion of 
agents in state 9 = A as they can bind to the membrane. hi{r) and h2{r) determine the agent's 
motion along the microtubule filaments, ho^r) on the other hand represents the cell topology, 
i.e. it generates a repelling force close to the cell membrane and to the nucleus, but is is simply 
a constant inside the cell, because free diffusion inside the cell should not be affected. 

The only time-dependent component of the effective field is /i3(Ar,t), which affects the fusion 
processes between vesicles (fission neglected). In fact, this is a short range attraction potential 
which increases with decreasing distance Ar = \ri — rj\. Since agents move, /13 changes in time 
depending on their actual positions ri{t) and internal states 9i(t). 

In order to describe how the effective potential results from the different field components, we 
have to consider the response function a[6i{t)] that determines which of the field components are 
actually "seen" by the agents conditional on their internal states. In accordance with the above 
distinction, we specify: 

a[ei{t)\ h''{r,6,t) = aoho{r) + Se„i aihi{r) + 6e„2 a2h2{r) 

+ <5e,,4 a^hAir) + 6g^^3 a3/i3(Ar, t) (4) 

The different ak are dimensional constants. With eqn. Q the motion of every agent is specified 
according to eqn. ([T]). We point out that agents with 9i = behave like simple Brownian particles 
which are not affected by any conditions inside the cell, except for the boundary conditions. 

Master equation The state of each individual agent is now described by a triple of three 
different state variables {ri{t),9i{t), Si{t)} that can change in time according to the processes 
specified above. The multi-agent system is thus described by the grand- canonical A^-particle 
distribution function 

Pn = PNir,9,s,t) = PN{ri,9i,si,...,rN,9N,SN,t), (5) 
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which describes the probability density of finding N Brownian agents with the distribution of 
internal parameters 9, positions r and sizes s at time t. Note that N is not constant but can 
change over time due to fusion (and fission) events. 

The complete dynamics for the ensemble of agents can be formulated in terms of a multivariate 
master equation: 

N 



dt 



PN{r,9,s,t)=-Y,\v 



i=l 
N 



a[ei{t)] 



70 



V,h%r,t,9)PMir,e,s,t) 



+ E E H(^^\(^'i) PN{r,9',,9\s,t)-u{9',\9,)PN{r,9,s,t)] 

i=i e[^e, 

N 

+ [w(si + Sj\si,Sj)PN+i{r^ ,9i = 3,9j = 3,9*,s*,t) 

i=l i<j 

-u}{si,Sj\si + Sj)PN{r,9i = 3,9j = 3,9, s,t)] 



DAiPN{r,9,s,t)} (6a) 
(6b) 



(6c) 



The first part of the multivariate master eqn., (6a), describes changes of the probability distribu- 



tion due to movements of agents either by diffusion or following some gradients of the effective 



field. The second part, (6b), considers all possible changes of the distribution of internal states. 



9, where 9* denotes "neighboring" states that differ from 9 only by the element explicitly given. 



The third part, (6c), eventually describes the fusion process by any two agents. This leads to a 
"gain" if the total number of vesicles is decreased hy N + 1 ^ N and distribution changes result 
in r* — )• Zli s* — )• s, or to a "loss" for any other process with N ^ N — 1. 

The multivariate master equation has the advantage of considering all possible processes on the 
agent level in a stochastic framework. After completely specifying the transition rates involved. 



we are able to solve this equation by means of stochastic computer simulations (see Section 4.4). 

We note that from eqn. ([6|, one can in principle derive a macroscopic density equation by 
introducing the agent density of the grand-canonical ensemble: 



n{r,t)= dri...drN~iPN{ri, ...,rN-i,r,t) 

N=l '' 



(7) 



The resulting equation would have the known structure of a reaction-diffusion equation for n(r, t). 
While this may be the "classical" way of investigation, we point out that for the agent-based 
approach proposed here the stochastic framework is the more appropriate one because it refers 
to the individual processes on the agent level. We emphasize that the agent-based approach is 
better suited for computer simulations of spatiotemporal patterns because it provides a stable and 
fast numerical algorithm. This becomes especially important in the case of large density gradients, 
which may considerably decrease the time step allowed to integrate the related dynamics. 
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4 Setup for stochastic simulations 

4.1 Velocities and transition rates 

According to the distinction between the activity modes given in Table [T| the movement of 
the agents can be of two types: free motion by diffusion processes, described by the diffusion 
coefficient D, and hound motion along the microtubuh. Both types are described by eqn. ([T]). In 
order to calculate the different Vi{6), we would need to specify explicitly the related components 
of the effective potential, he{r,t) that result in the directed motion along the microtubuli. To 
simplify the procedure and match it with experimental findings, we may instead consider the 
value of Vi{9) as given by experiments. Then, the role of the field component is reduced to simply 
keep the agents on the microtubuli as long as they are in the respective internal states 6. I.e. 
both Vi/ii(r) and Vj/i2(T') just determine the direction of motion along the microtubulus to 
which agent i is attached, whereas /lo(^) specifies the boundary condition given by the location 
of the cell nucleus and the cell membrane (see also Figure |2]). Table [2] provides the values for the 
velocities which are assumed to be constant and the same for all agents in the respective internal 
state. 



State 


Parameter value 


Description 


^ = 


D = 10"^^ mVs 


Diffusion coefficient in cytoplasma 


9 = 1 


vi = 0.33 /xm/s 


Velocity of vesicle on MT towards plus-ends 


6 = 2 


V2 = 0.33 //m/s 


Velocity of vesicle on MT towards minus-ends 



Table 2: Parameters describing the free {6 = 0) or bound {6 = 1,2) motion of agents. MT 
means 'microtubulus'. Note that the diffusion coefficient is related to the friction coefficient via 
7o = ksT/D, which results in 70 = 4 • 10~^kg/s. Values according to j5]. 

As a next step, we need to specify the transition rates uj{0'\0) between different modes of activ- 
ity. Again, instead of providing explicit expressions, we may simply take values obtained from 
experiments. These values are available to us only for the unperturbed state of the cell, i.e. for a 
baseline or reference case in which conditions inside the cell are not changed on purpose. In the 
following, baseline values are indicated by the superscript h. Table [3] lists all biologically relevant 
transition rates for the unperturbed scenario, as obtained either from the literature or from own 
experiments. 

4.2 Reduced transition rates 

The given set of transition rates leaves us with a large degree of freedom which however turns out 
to be a pitfall: if we wanted to compare the simulations with patterns observed in the experiment, 
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Parameter Value 


Description 




0)=0.05 


Transition from diffusive to MT (plus-end) bound state 




1)=0.33 


Transition from MT-bound (plus-end) to diffusive state 


u:\2 


0)=0.05 


Transition from diffusive to MT (minus-end) bound state 


a;^(0 


2)=0.33 


Transition from MT (minus-end) bound to diffusive state 




1)=0.01 * 


Transition from MT-bound to MT-bound state with fusion 


u:\l 


3)=0.02 " 


Transition from MT-bound state with fusion to MT-binding 




2)=0.01 * 


Transition from MT-bound to MT-bound state with fusion 




3)=0.02 * 


Transition from MT-bound state with fusion to MT-binding 




4)^0.0025-0.0033 s'^ 


Internalization rate of vesicles 




0)^0.00083-0.0012 s^^ 


Recycling rate of vesicles 



Table 3: Biologically relevant transition rates iJ'{6'\9) for the baseline case (unperturbed con- 
centrations of proteins). Values according to [5], values with * are from own experiments, see 

[IJ. 



we would need to raster the full parameter space to find appropriate combinations of transition 
rates which lead to a realistic outcome. There are basically two ways to reduce this parameter 
space: (i) we use known values of the transition rates as e.g. reported in the literature, (ii) we 
introduce reduced transition rates, assuming that not all processes are really independent. We 
follow a combination of the two, defining the reduced transition rates as given in Table |4] 



9.1 = 


U}{1 


0) 


Affinity for microbule plus-ends 


LO{0 


1) 


^2 = 




0) 


Affinity for microbule minus-ends 


io{0 


2) 


^3 = 


u{3 


1) 


Fusion tendency on microbule plus-ends 


oo{l 


3) 


= 


a;(3 


2) 


Fusion tendency on microbule minus-ends 


u}{2 


3) 


^5 = 




4) 


Internalization versus recycling rate 


a;(4 


0) 



Table 4: Reduced transition rates fij, to combine two transition rates u which refer to the 
same intracellular (transport) mechanism but describe inverse processes. 

In the following we assume that = ^^4^ = ftp which denotes the fusion tendency that is 
supposed to be independent of the kind of motor protein involved (this was kinesin for 0,^ and 
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dynein for $74). With this, we finally have reduced the transition rates to 4 parameters given by 
r^i, il2, and r^s. 

These reduced parameters of course do not fully determine the individual transition rates which 
are needed to recover the correct dynamics. Therefore, for one of the transition rates we use the 
baseline value from the literature, as given by Table [s] So, we define e.g. Q.i = a;(l|0)/a;^(0|l) or 
VLi =uj\l\Q)/uj{^\l). 

Eventually, we want to emphasize the distinction between the reduced transtion rates fi" for 
the unperturbed case (control cell) and for the perturbed case, where the protein k was 
manipulated. Similar to the discussion in Section [4.1[ they refer to the same transition rates w, 
but with different concentrations Cfc. 

4.3 Boundary and initial conditions 

In order to complete the setup for the stochastic computer simulations, we still need to specify 
the boundary conditions which refer to the cell geometry. Figure [2] (a,b) shows the two different 
cell geometries chosen, a rather circular cell and an elongated one. The outer boundary of the 
cell membrane and the inner boundary of the nucleus are both assumed to be impermeable walls, 
described by a hard sphere potential hQ{r). 

The interior of the cell contains the cytoskeleton, i.e. both actin filaments and microtubuli (MT) 
which provide boundary conditions for the motion of vesicles (see also Figure [T| . As already 
discussed, we do not consider motion along actin filaments because the related processes do not 
contribute to vesicle pattern formation. MT, on the other hand, are abundant and always point 
to the microtubule organization center (MTOC) which is assumed to be in the perinuclear area 
on the side with the largest part of the cytosol (see Figure |2] a, b). 

Regarding the initial conditions, we need to specify the number, the position, the internal state 
and the size of the agents at t = 0. For our model, we assume that initially all vesicles are bound 
to the cell membrane, i.e. at t = agents start with 0j(O) = 4 at a position rj(0) randomly 
chosen from the cell boundary. In order to determine the initial number of agents, we start 
from the experimental observation that an average cell (of the type considered) contains about 
200 internalized vesicles at steady state. This number excludes vesicles still bound to the cell 
membrane. We further know from experiments that internalization rates, i.e., transition rates 
from the initial bound into the free moving state (and vice versa), i.e. a;(0|4) = 3.3 • 10~'^s~^ and 
(^(4|0) = 8.3 • 10~^s~"^ (see Table [3|. If A^(0) = A'4(0) denotes the (unknown) number of agents 
at t = (all assumed to be in the bound state) and A^(0) — A'l* = 200 the (known) number of 
agents no longer in the bound state at steady state (st : steady state), then we can postulate for 
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Figure 2: (left) (a) Circular geometry of a non-treated (control) cell; (right) (h) Elongated ge- 
ometry of a perturbed cell (treated with siRNA, which leads to the depletion of CDK8). The 
scale bar corresponds to h(j.m. The microtubule filaments are schematically sketched, pointing 
toward the MTOC (microtubule organization center). Vesicle in an internal state 6 G {1,2,3} 
move along microtubules in both directions. 



the dynamics for N4^{t) the following rate equation: 

= -^(0|4) N^) + ^(4|0) [iV4(0) - N^)] (8) 
After a time t ~ l/[a;(0|4) -|-6t;(4|0)], this dynamics reaches a steady-state solution 

from which we can calculate A'^4(0) assuming that N{{))—N^ = 200. This approximation neglects 
the fact that the total number of vesicles have decreased at time t because of the fusion between 
vesicles. Thus, we may slightly increase the initial number of agents and have eventually chosen 
iV4(0) = 350. 

It remains to determine the initial distribution of vesicle sizes. We want to start from a most 
realistic one because (a) later we want to compare the time scale of structure formation with the 
experimental observation, and (b) because our modeling setup has neglected the fragmentation 
rate of vesicles (which would be needed if an arbitrary initial distribution needs to relax into 
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a realistic one). Again, we rely on experimental observations [T] that have found a log-normal 
distribution of vesicle sizes: 

where /i and a are the mean and standard deviation of the variable's natural logarithm. 
4.4 Stochastic Simulation Technique 

We now have determined all ingredients for stochastic computer simulations which include the 
following dynamic processes (specified on the agent level): 

Initialization At t = to, 350 agents with 6i{0) = 4 are randomly placed at the cell boundary. 



Their initial size Sj(0) is drawn from the log-normal distribution, (10). 



Movement Agents can change from the bound state into the free moving state at a rate Ci;(0|4) 
and from the free moving state into directed motion at rates a;(l|0), a;(2|0). They all move 
according to the equation of motion, ([T]). 

In our modeling approach, only agents with the internal states 9 £ {1,2,3} move along 
the MT. We assume that, whenever an agent switches from 9 = into either 9 = 1 or 
9 = 2, i.e. from a free motion into a bound motion, a MT is "always at hand". If the agent 
switches from 9 S {1,2} into 9 = 3 where it is ready to fuse, it continues to move into the 
same direction as before, until it collides with another agent. 

Fusion Precisely, fusion occurs only on MT. Agents need to be in state 9 = 3 and in a sufficiently 
close distance. After fusion, the smaller agent "disappears", whereas the larger one has 
increased its size, but keeps the internal state ^ = 3 until one of the transitions a;(2|3) or 
w(l|3) happen. 

Inversion Agents bound to the MT can become freely moving at the rates a;(0|2), a;(0|l), 
whereas free moving agents can be bound again to the cell membrane at the rate Ct;(4|0). 

Concurrency We assume that agents can move and transit into different internal states at the 
same time. This allows us to decouple the motion of agents from the various reactions 
(change of internal states and fusion). 

The state of the multi-agent system is at any time completely described by the A^-particle 
distribution function Pj\f (r, 9, s, t) . However, because of the dynamical processes, the system 
state always changes and its average "life time" Tm is just the inverse of the sum over all possible 
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transition rates that can change the given state (including changes of position, internal states, 
and sizes). 

Because we have to deal with movement and reactions at the same time, we have chosen a 
sufficiently small fixed time interval At = 0.02s to solve the equations of motion of each individual 
agent. To answer the question if in the respective time interval also a change of the internal states 
or a fusion process occurs, we proceed as follows: 

Each of the possible reactions has a different probability to occur, which is determined by the 
respective transition rate and the available time, At. In order to pick one of the possible reactions, 
we draw a uniformly distributed random number U G RND[0, 1] and choose the process z that 
satisfies the condition 

2 2 + 1 

^w,(-)At<C/<^L^,(-)At (11) 

Because At was chosen such that Ylm=i^n{')^'t < 1, none of the possible processes is excluded 
from being picked. On the other hand, it may occur that none of the processes is being chosen if 
the sum is much smaller than 1 and U close to 1. Then no reaction occurs during the respective 
time interval, but movements take place. 

5 Results 

5.1 Spatio-temporal vesicle patterns 

Figure [3] presents computer simulations of the vesicle patters for both the unperturbed (control) 
cell and the perturbed cell (cf also Figure |2]) . We emphasise that these simulations refer to real 
spatial and time scales, so they should be comparable, at least in a statistical sense, to patterns 
observed from experiments. These experiments are reported in [1] and have motivated the choice 
of the reduced transition rates, which are treated in this paper as free parameters. 

Comparing the pattern formation in the perturbed cell with the one in the control cell, we note a 
number of differences: we observe a localization of large vesicles on the one hand at the tips of the 
elongated branched-out perturbed cell, on the other hand large vesicles are as well localized in 
the perinuclear area of this cell. In the unperturbed cell, the majority of large vesicles are located 
around the nucleus and vesicles are spread over the entire cell surface getting more sparse towards 
the cell periphery. 

Comparing the vesicle size distribution of both the perturbed and the unperturbed cell, we find 
that they preserve the form of the log-normal distribution, but the mean value |U, in the course 
of time, shifts to significant larger values in the perturbed cell (see also Table |5]). The perturbed 
cell displays a geometry that may have facilitated fusion of vesicles in its branches within which 
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they accumulate. The transport of vesicles from these branches to the nucleus of the perturbed 
cell seems to be suppressed. 



t 


1 min 


5 min 


10 min 


15 min 




2.69 


2.82 


2.94 


2.96 


cr" 


0.89 


0.94 


0.94 


0.96 




2.87 


2.99 


3.00 


3.04 


a 


0.90 


0.90 


0.94 


0.86 



Table 5: Mean and standard deviation a{t) of the log-normal distribution at different 
times t (min) for the unperturbed (superscript u) and the perturbed cell. 



5.2 Estimating concentration dependence 

As pointed out above, we found very realistic vesicle patterns both for the perturbed and the 
control cell for the given set of reduced transition rates (see Figure [3| . These transition rates 
are treated as free parameters in our model - but provided they are correct (which can only 
be confirmed by comparing the patterns with experiments, statistically) they allow an indirect 
estimation of the concentration dependence of the transition rates oj[9'\9). 

As already stated, the transition rates are given to us only for the unperturbed state of the cell. 
Table |3]presents the values for the baseline case. One should note, however, that the baseline value 
not necessarily describes the experimental situation because it was obtained under conditions 
which are hardly reproducible completely. If we for example change the concentration Ck of some 
protein k inside the cell which is involved in processes of fusion or directed motion, this will 
certainly change the value of the respective transition rates, i.e. oj = a;(cfc, q, ...). Hence, it is not 
only sufficient to know the values of the baseline case, we should also know how these values of 
the transition rates change with the concentrations c^. If we denote the transition rates in the 
unperturbed case by lj" (omitting the 6 dependence at the moment) and the respective protein 
concentrations by c" we may assume in first-order approximation the following expansion: 

9a; 



^{ck, ci^k = cf^k) = w"(cfc, cf^k) + 



ick-4). (12) 



To further specify the functional dependency uj{ck) we make the following ansatz: 

(^{ck,ci^k = cf^k) = ■ exp (^Kk '^^ ^ . (13) 

which satisfies (jj{ck,ci^k = c^fc) = '^^ for Ck = c^- The important parameter Kk denotes the 
impact that a change of concentration Ck has on the respective transition rate, i.e. it is a measure 
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1 min 5 min 10 min 15 min 




Figure 3: Computer simulations of vesicle patterns on real time scale (in min) and spatial scale 
(scale bar: 5^m ). (top) unperturbed cell: {Qi, f22, ^2^?, Q5}={0. 925, 1.1, 1.0, 0.85}, (bottom) 
perturbed cell: {^i, Vt2, ^^5}={1-3, 1.2, 1.2, 0.55}. The histograms show the evolution of the 
vesicle size distribution together with the fitted log-normal distribution (red line). Values for /i 
and <T are given in Table |5j 

of sensitivity toward that particular protein. Of course, = K,k{0'\0) in full notation, i.e. the 
value does not only change across proteins, but the impact also changes for different transitions. 
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Putting eqs. (12), (13) together, we arrive at: 



Awfc = w(cfc,Q^fc = cf^fc) - = t^"^(cfc - cl) (14) 

In many experimental cases, as e.g. in RNAi screens, the perturbation of a protein concentration 
leads to ct — ?• 0, i.e. we are interested in the dynamics in the absence of a given protein. With 
this assumption we finally have: 

^0Jk{e\d') = -Kk{e\0') a;"(0|0')- (15) 

This allows us to relate two different dynamical scenarios and their respective outcome in terms 
of the vesicle patterns: (i) the unperturbed scenario with experimentally known concentrations 
and known transition rates, and (ii) the perturbed scenario^ where different proteins may be 
absent. 

As an example, let us investigate how changes in the concentration of the protein k = CDK8 
affect the transition rates a;cDK8(l|0) ^-nd wcDK8(0|l)- Given the parameters in Figure |3] the 
reduced transition rates return fif — = 0.375, where refers to the case where the 

protein concentration Cfc = ccdks = 0) whereas fi" refer to the unperturbed case. Dividing eqn. 



(15) by a;''(0|l), we find 



Af^f°^« = - 0^ = 0.375 = • KCDK8(1|0) 

= -^^1 • KCDK8(1|0), (16) 

from where it follows that kcdK8(1|0) ~ —0.41, which describes the sensitivity toward changes 
in the concentration of CDK8. 

Knowing the difference Afif^^^, the transition rates Uk are not fully determined because of 
= w"(l|0)/a;*(0|l), 0^°^*^ = wcdK8(1|0)/w''(0|1). Hence, we can now discuss two different 
cases which refer to two hypotheses about the transport of vesicles towards microtubule minus- 
ends. 

The first hypothesis of our model states that the transition rate a;cDK8(l|0) in cells that are 
silenced for CDK8 reads: 

^cdk8(1|0) = .."(110) • exp f - 0.41 • ^cdk8_^^\ _ ^^^^ 

Assuming ccdks ~ 0, it follows, that ti;cDK8(l|0) is increased by a factor of approximately 1.5 
with respect to a;"(l|0). This means that silencing the protein CDK8 increases the transition of 
vesicles freely diffusing in the cytosol to vesicles being transported towards the MT-plusends. We 
thus hypothesize that CDK8 is either directly or indirectly involved in the docking of vesicles on 
microtubule filaments via kinesins. 
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On the other hand, kcdk8(0|1) can be obtained by the following relation: 




= fi^^^«-KCDK8(0|l), 



(18) 



from where it follows that AvcdK8(0|1) ~ 0.29. The second hypothesis of our model states that 
the transition rate a;cDK8(0|l) in cells that are silenced for CDK8 reads: 



Setting CCDK8 ~ Oj we find that a;cDK8(0|l) is decreased by a factor of approximately 0.75 with 
respect to u;"(0|l). This means that silencing the protein CDK8 leads to a decrease in transition 
of vesicles transported towards the MT-plusends to vesicles freely diffusing vesicles. We thus 
hypothesize that CDK8 may as well be involved in the undocking of vesicles bound to MT 
filaments via kinesins. 

In a similar manner, we could estimate the concentration dependence of other transition rates, 
using the scaled transition rates {9!^^^^ ,0.^^^^) and (fif , O^, fi^). This enables us to 
further develop a number of hypotheses regarding the effect of silencing the protein CDK8 on 
the intracellular transport. 

6 Discussion 

6.1 Motivation of the model 

Genomic and pharmaceutical research nowadays heavily relies on systematic screens in which the 
perturbation or silencing of specific proteins affects the abundance of vesicle patterns observed 
in a population of cells. The study of vesicle pattern formation thus is important to improve 
our understanding of the function of genes. To learn how vesicle patterns are formed, we have 
set up an agent-based model of intracellular transport inside a single cell. Agents represent 
vesicles which move either by diffusion in the cytosol or are transported along the cytoskeletal 
filaments through molecular motors. Vesicles further interact with other vesicles by fusion or 
fission, or they interact with the cytoskeletal filaments, the cell membrane and the nucleus. This 
interaction is controlled and regulated by specific proteins which are synthesised inside the cell by 
transcriptionally active genes. The activity of these genes thus represents the control parameters 
of our system. 

Treating vesicles as Brownian agents with an internal degree of freedom allows to formally 
derive a model that captures all relevant processes in the formation of vesicle patterns. Five 



WCDK8(0|1) =w"(0|l)-exp 0.29- 



CCDK8 - CCDK8 
^CDK8 



) 



(19) 
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different values of the internal degree of freedom define the vesicle's different modes of activity. 
Transitions between these modes occur as stochastic events related to the binding of proteins 
to the vesicle's coat, or to signalling events. Proteins involved in those processes can control 
the vesicle's activity which in turn determines the process of pattern formation. We therefore 
assume that the transition rates between different value of the internal degree of freedom depend 
on the concentration of a cell's proteins. To determine the precise transition rates as functions of 
protein concentrations would require the full knowledge about the regulatory structure of gene 
networks. We simplified this situation by assuming that the transition rates can be decomposed 
into a product of the unperturbed transition rates and an exponential function depending on 
the concentration of single proteins. We further simplified our model by assuming that the cell 
membrane, the cytoskeletal filaments and the nucleus are static and that the diffusion coefficient 
and velocities along cytoskeletal filaments are constant parameters of the model. In contrast, 
the transition rates represent the free parameters of our model and have been varied. We have 
reduced the number of 10 original transition rates to 4 scaled transition rates. This introduced 
an ambiguity in the interpretation of the simulation results with respect to the original transition 
rates. We could cope with this situation by deriving two complementary hypothesis about the 
influence of the proteins involved, which could be tested experimentally. 

6.2 Possible comparison with experiments 

In order to relate our computer simulations to reality, experimental data are needed to calibrate 
the simulations. Whenever available, baseline values for the transition rates obtained from exper- 
iments have been included. This allowed us to generate patterns on real time and spatial scales. 
From every simulated pattern, four features can be obtained: size, relative distance to nucleus, 
number of vesicles within a fixed radius around each vesicle and number of vesicles per cell area. 
To measure the dissimilarity between the simulated pattern and the experimentally observed one, 
it is useful to compute the symmetrized KuUback-Leibler divergence of the two corresponding 
vesicle feature distributions [T] , to find out for which parameters the simulated pattern provides a 
minimum divergence to the experimentally observed vesicle patterns. An interpretation of these 
findings in comparison with hypotheses generated from the computer simulations allows to draw 
conclusions about the underlying processes, in particular about the role of the genes involved. 

A comparison of simulated and experimentally measured vesicle patterns also involves a dimen- 
sional problem: our simulations are performed in 2d, whereas experimental patterns result from 
vesicle motion in 3d. In principle, this would require to correct for the distances as consequence 
of the projection from 3d to 2d. Since mammalian cells are relatively flat with the exception of 
the nucleus area, we have omitted this correction. But we considered the fact that vesicles could 
pass below/above each other by not assuming mutual exclusion in our 2d simulations. 
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6.3 Future directions 



We emphasize again that our computer simulations lead to testable hypotheses about the influ- 
ence of specific genes, such as CDK8, on the formation of vesicle patterns. One future apphcation 
of our model of intracellular transport concerns the silencing of multiple genes. Let us denote 



with i and / two genes which are silenced, then according to ansatz in eqn. (13) we find for the 
transition rates 

u:r0,\e,) = • exp I ^,,{e'^\e + ^0^\e^f-i^\ . (20) 

Once we have determined Ki{9'-\9j) and ni[6'-\6j) from single gene silencing experiments and their 
related simulations of our model of intracellular transport, we can predict the resulting pattern 



and dynamics of the combined silencing of these two genes. However, (eq. 20) is only valid, if 
genes i and / are not in the same (regulatory) pathway. Such a prediction represents an in-silico 
experiment and can be tested experimentally. 



Acknowledgment 

M.B. could benefit from numerous stimulating discussions with Lucas Pelkmans. 



References 

[1] Birbaumer, M. (2010). Statistical analysis and modeling of endocytic vesicle pattern forma- 
tion. Ph.D. thesis, ETH Zurich, Disseration Nr. 19322. 

[2] Conner, S.; Schmid, S. (2003). Regulated portals of entry into the cell. Nature 422(5875), 
37-47. 

[3] Dautry-Varsat, A. (1986). Receptor-mediated endocytosis: the intracellular journey of trans- 
ferrin and its receptor. Biochimie 68 , 375-81. 

[4] Dinh, A.-T.; Pangarkar, C; Theofanous, T.; Mitragotri, S. (2006). Theory of Spatial Pat- 
terns of Intracellular Organelles. Biophysical Journal 90(10), L67-L69. 

[5] Dinh, A.-T.; Pangarkar, C; Theofanous, T.; Mitragotri, S. (2007). Understanding Intracel- 
lular Transport Processes Pertinent to Synthetic Gene Delivery via Stochastic Simulations 
and Sensitivity Analyses. Biophysical Journal 92(3), 831-846. 

[6] Doherty, G.; McMahon, H. (2009). Mechanisms of Endocytosis. Annual Review of Biochem- 
istry 78, 857-902. 



20/: 



21 



Mirko Birbaumer, Frank Schweitzer: 
Agent-Based Modeling of Intracellular Transport 
European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 



[7] Ebeling, W.; Schweitzer, F. (2003). Self-Organisation, Active Brownian Dynamics, and 
Biological Applications. Nova Acta Leopoldma NF 88(332), 169-188. 

[8] Ebeling, W.; Schweitzer, F.; Tilch, B. (1999). Active Brownian particles with energy depots 
modeling animal mobility. Biosystems 49(1), 17 - 29. 

[9] Fiore, P. D.; Camilli, P. D. (2001). Endocytosis and Signaling: An Inseparable Partnership. 
Cell 106(1), 1 - 4. 

[10] Garcia, V.; Birbaumer, M.; Schweitzer, F. (2011). Testing an agent-based model of bacte- 
rial cell motility: How nutrient concentration affects speed distribution. European Physical 
Journal B 82(3-4), 235-244. 

[11] Hannon, G. (2002). RNA interference. Nature 418(6894), 244-251. 

[12] Mach, R.; Schweitzer, F. (2007). Modeling Vortex Swarming In Daphnia. Bulletin of Math- 
ematical Biology 69, 539-562. 

[13] Marsh, M.; Helenius, A. (2006). Virus entry: open sesame. Cell 124(4), 729-740. 

[14] Mercer, J.; Helenius, A. (2008). Vaccinia Virus Uses Macropinocytosis and Apoptotic 
Mimicry to Enter Host Cells. Science 320(5875), 531-535. 

[15] Roth, M. (2006). Clathrin-mediated endocytosis before fluorescent proteins. Nat Rev Mol 
Cell Biol 7, 63-68. 

[16] Schweitzer, F. (2003). Brownian Agents and Active Particles: Collective Dynamics in the 
Natural and Social Sciences. Springer, Berlin. 

[17] Schweitzer, F.; Garcia, D. (2010). An agent-based model of collective emotions in online 
communities. European Physical Journal B 77(4), 533-545. 



21/: 



21 



